# fig1-simulate-data.R
# 4/17/24
setwd("~/Dropbox/Benton-Jordan-Philips/final-JOP-replication")
library(rugarch)
library(readstata13)
library(MASS) # for mvnorm
library(dynamac)
library(reshape2)
library(ggplot2)
library(Matrix)
library(foreign)

# Variance. all of these are set up to have multiple values in the loop
# het equation
the.hetx.mean <- 0	# This is for het_x
the.hetx.sd <- c(0.5, 1)		# Again, X and Z comparable
the.beta.hetx <- c(0.1, 1)    # Make X and Z have comparable effects

the.hetz.mean <- 0
the.hetz.sd <- c(0.5, 1)
the.beta.hetz <- c(0.1, 1)

# Limit arch/garch/beta.hetx to good range
the.arch <- c(0.2)
the.garch <- c(0.2, 0.5, 0.8)

# Simulation parameters. Seeds and lengths.
# Number of sims
num.seeds <- 10

set.seed(43)
the.seeds <- round(runif(num.seeds, 1, 1000000))
unique(the.seeds)
the.lengths <- c(400) # actually going to add 100 later for burnin in loop

# burnin
burnin <- 100

# Common mean parameters + x_mean. None of these are set up to have multiple values in the loop.
beta.x <- 2 # coef on x in mean equation
beta.lx <- -1 # coef on l.x in mean equation
phi <- 0.1 # coef on ldv in mean equation
const <- 2 # constant needed to have a value of mean.y
omega <- 0.5
# And one shared series for the x_mean
set.seed(1000)
x.mean.full <- rnorm(max(the.lengths) + burnin)

# Columns are the number of X*2 (creating a series of x and y), shared xmean, shared beta.hetz, shared time, a column for combo, and the 
#  combo parameters (arch, garch, beta.hetx, seeds, lengths, means, and sds)
bigmat <- matrix(rep(NA, 
	(max(c(max(the.lengths), (length(the.arch)*
								length(the.garch)*
								length(the.beta.hetx)*
								length(the.hetx.mean)*
								length(the.hetx.sd)*
								length(the.beta.hetz)*
								length(the.hetz.mean)*
								length(the.hetz.sd)*
								length(the.seeds)*
								length(the.lengths))))) * # this is the rows
		(((length(the.arch)*
								length(the.garch)*
								length(the.beta.hetx)*
								length(the.hetx.mean)*
								length(the.hetx.sd)*
								length(the.beta.hetz)*
								length(the.hetz.mean)*
								length(the.hetz.sd)*
								length(the.seeds)*
								length(the.lengths)) * 3) + 13)), # this is the columns. *3 for x/y/z, 13 for shared
	nrow = (max(c(max(the.lengths), (length(the.arch)*
								length(the.garch)*
								length(the.beta.hetx)*
								length(the.hetx.mean)*
								length(the.hetx.sd)*
								length(the.beta.hetz)*
								length(the.hetz.mean)*
								length(the.hetz.sd)*
								length(the.seeds)*
								length(the.lengths))))))

# Fill the first column after X/Y with time
bigmat[(1:max(the.lengths)),(((length(the.arch)*length(the.garch)*
						length(the.beta.hetx)*length(the.hetx.mean)*length(the.hetx.sd)*
						length(the.beta.hetz)*length(the.hetz.mean)*length(the.hetz.sd)*
						length(the.seeds)*length(the.lengths)) * 3) + 1)] <- seq(1, max(the.lengths), 1)
# Fill the second column after X/Y with shared x_mean series
bigmat[(1:max(the.lengths)),(((length(the.arch)*length(the.garch)*
						length(the.beta.hetx)*length(the.hetx.mean)*length(the.hetx.sd)*
						length(the.beta.hetz)*length(the.hetz.mean)*length(the.hetz.sd)*
						length(the.seeds)*length(the.lengths)) * 3) + 2)] <- x.mean.full[(burnin + 1):(max(the.lengths) + burnin)]					

# Incrementers to fill matrix
combo <- 1

create.data <- proc.time()

for(u in 1:length(the.seeds)) {
	for(v in 1:length(the.hetz.mean)) {
		for(w in 1:length(the.hetz.sd)) {
			for(x in 1:length(the.hetx.mean)) {
				for(y in 1:length(the.hetx.sd)) {
					set.seed(the.seeds[u])
					z.het.full <- rnorm((max(the.lengths) + burnin), mean = the.hetz.mean[v], sd = the.hetz.sd[w])
					x.het.full <- rnorm((max(the.lengths) + burnin), mean = the.hetx.mean[x], sd = the.hetx.sd[y])
					y.full <- error <- sigma2 <- rep(NA, length(x.mean.full))
					y.full[1] <- const/1-phi
					sigma2[1] <- error[1] <- rnorm(1)
					for(m in 1:length(the.arch)) {
						for(n in 1:length(the.garch)) {
							for(o in 1:length(the.beta.hetx)) {
								for(p in 1:length(the.beta.hetz)) {
									for(i in 2:length(y.full)) {
										sigma2[i] <- exp(omega + the.beta.hetx[o]*x.het.full[i] + the.beta.hetz[p]*z.het.full[i]) + the.arch[m]*(error[i - 1]*error[i - 1]) + the.garch[n]*sigma2[i - 1]
										error[i] <- rnorm(1, 0, sqrt(sigma2[i]))
		  								y.full[i] <- const + phi*y.full[i - 1] + beta.x*x.mean.full[i] + beta.lx*x.mean.full[i - 1] + error[i]
									}
									for(z in 1:length(the.lengths)) {
										# Fill X (leave gaps for all x), then fill y. disregard burnin
										bigmat[(1:the.lengths[z]), combo] <- x.het.full[(burnin + 1):(burnin + the.lengths[z])]
										
										bigmat[(1:the.lengths[z]), (((length(the.arch)*length(the.garch)*
											length(the.beta.hetx)*length(the.hetx.mean)*length(the.hetx.sd)*
											length(the.beta.hetz)*length(the.hetz.mean)*length(the.hetz.sd)*
											length(the.seeds)*length(the.lengths)) * 1) + combo)] <- y.full[(burnin + 1):(burnin + the.lengths[z])]
											
										bigmat[(1:the.lengths[z]), (((length(the.arch)*length(the.garch)*
											length(the.beta.hetx)*length(the.hetx.mean)*length(the.hetx.sd)*
											length(the.beta.hetz)*length(the.hetz.mean)*length(the.hetz.sd)*
											length(the.seeds)*length(the.lengths)) * 2) + combo)] <- z.het.full[(burnin + 1):(burnin + the.lengths[z])]
											
										# columns 1 and 2 and 3 after the sims are time and shared x.mean 
										# Column 3. the.beta.hetx		
										bigmat[combo, (((length(the.arch)*length(the.garch)*
											length(the.beta.hetx)*length(the.hetx.mean)*length(the.hetx.sd)*
											length(the.beta.hetz)*length(the.hetz.mean)*length(the.hetz.sd)*
											length(the.seeds)*length(the.lengths)) * 3) + 3)] <- the.beta.hetx[o]
										# Column 4. the.hetx.mean
										bigmat[combo, (((length(the.arch)*length(the.garch)*
											length(the.beta.hetx)*length(the.hetx.mean)*length(the.hetx.sd)*
											length(the.beta.hetz)*length(the.hetz.mean)*length(the.hetz.sd)*
											length(the.seeds)*length(the.lengths)) * 3) + 4)] <- the.hetx.mean[x]
										# Column 5. the.hetx.sd
										bigmat[combo, (((length(the.arch)*length(the.garch)*
											length(the.beta.hetx)*length(the.hetx.mean)*length(the.hetx.sd)*
											length(the.beta.hetz)*length(the.hetz.mean)*length(the.hetz.sd)*
											length(the.seeds)*length(the.lengths)) * 3) + 5)] <- the.hetx.sd[y]						
										# Column 6. the.beta.hetz		
										bigmat[combo, (((length(the.arch)*length(the.garch)*
											length(the.beta.hetx)*length(the.hetx.mean)*length(the.hetx.sd)*
											length(the.beta.hetz)*length(the.hetz.mean)*length(the.hetz.sd)*
											length(the.seeds)*length(the.lengths)) * 3) + 6)] <- the.beta.hetz[p]
										# Column 7. the.hetz.mean
										bigmat[combo, (((length(the.arch)*length(the.garch)*
											length(the.beta.hetx)*length(the.hetx.mean)*length(the.hetx.sd)*
											length(the.beta.hetz)*length(the.hetz.mean)*length(the.hetz.sd)*
											length(the.seeds)*length(the.lengths)) * 3) + 7)] <- the.hetz.mean[v]
										# Column 8. the.hetz.sd
										bigmat[combo, (((length(the.arch)*length(the.garch)*
											length(the.beta.hetx)*length(the.hetx.mean)*length(the.hetx.sd)*
											length(the.beta.hetz)*length(the.hetz.mean)*length(the.hetz.sd)*
											length(the.seeds)*length(the.lengths)) * 3) + 8)] <- the.hetz.sd[w]
										# Column 9. combo
										bigmat[combo, (((length(the.arch)*length(the.garch)*
											length(the.beta.hetx)*length(the.hetx.mean)*length(the.hetx.sd)*
											length(the.beta.hetz)*length(the.hetz.mean)*length(the.hetz.sd)*
											length(the.seeds)*length(the.lengths)) * 3) + 9)] <- combo
										# Column 10. arch
										bigmat[combo, (((length(the.arch)*length(the.garch)*
											length(the.beta.hetx)*length(the.hetx.mean)*length(the.hetx.sd)*
											length(the.beta.hetz)*length(the.hetz.mean)*length(the.hetz.sd)*
											length(the.seeds)*length(the.lengths)) * 3) + 10)] <- the.arch[m]
										# Column 11. garch
										bigmat[combo, (((length(the.arch)*length(the.garch)*
											length(the.beta.hetx)*length(the.hetx.mean)*length(the.hetx.sd)*
											length(the.beta.hetz)*length(the.hetz.mean)*length(the.hetz.sd)*
											length(the.seeds)*length(the.lengths)) * 3) + 11)] <- the.garch[n]
										# Column 12. seed
										bigmat[combo, (((length(the.arch)*length(the.garch)*
											length(the.beta.hetx)*length(the.hetx.mean)*length(the.hetx.sd)*
											length(the.beta.hetz)*length(the.hetz.mean)*length(the.hetz.sd)*
											length(the.seeds)*length(the.lengths)) * 3) + 12)] <- the.seeds[u]
										# Column 13. length
										bigmat[combo, (((length(the.arch)*length(the.garch)*
											length(the.beta.hetx)*length(the.hetx.mean)*length(the.hetx.sd)*
											length(the.beta.hetz)*length(the.hetz.mean)*length(the.hetz.sd)*
											length(the.seeds)*length(the.lengths)) * 3) + 13)] <- the.lengths[z]
										combo <- combo + 1		
										if(endsWith(as.character(combo), "000")) {
											print(combo)
										}
									}
								}
							}
						}
					}
				}
			}
		}
	}
}
# > proc.time() - create.data
#    user  system elapsed 
#   6.429   0.233   6.550 

colnames(bigmat) <- c(paste("x_het_", seq(1, (combo - 1), 1), sep = ""),
						paste("y_comb_", seq(1, (combo - 1), 1), sep = ""),
						paste("z_het_", seq(1, (combo - 1), 1), sep = ""),
						"time", "x_mean", "beta_hetx", "mean_of_hetx", "sd_of_hetx", 
							"beta_hetz", "mean_of_hetz", "sd_of_hetz", "combination", "arch", "garch", "seed", "length")

bigmatdata <- data.frame(bigmat)
head(bigmatdata)
tail(bigmatdata)

table(bigmatdata$beta_hetx, bigmatdata$beta_hetz, bigmatdata$garch, bigmatdata$seed)


write.dta(bigmatdata, "data/simulated_wide.dta")
